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SUMMARY 

This  final  technical  report  summarizes  the  research  carried  out 
under  sponsorship  of  AFOSR  contract  number  F49620-88-C-0060 , 
awarded  to  AEDAR  Corporation  for  the  period  May  1,  1988  -  July 
31,  1989.  The  main  thrust  of  the  research  has  been  the 
establishment  of  feasible  approaches  for  using  computer  algebra 
to  derive  the  Green's  Function  of  interconnected  structural 
systems.  The  accomplishments  of  the  research  effort  during  this 
period  include  the  formulation,  analysis  and  coding  of  computer- 
algebraic  techniques  for  the  derivation  of  the  Green's  functions 
of  structural  systems  consisting  of  uniform  Euler-Bernoulli  beams 
which  are  either  supported  elastically  or  to  which  discrete 
substructures  have  been  attached.  The  computer-algebraic 
environment  used  was  MACSYMA.  Using  a  general  representation  of 
the  Green's  function  of  elastic  structures,  the  process  of 
deriving  the  Green's  function  of  such  interconnected  structures 
consisted  of  the  followina  staqes:  (1)  Solution  of  the  integral 
equation  which  determines  the  transfer  function  of  the  updated 
system,  (2)  Determination  of  the  characteristic  parameters  of  the 
transfer  function  for  the  updated  system,  (3)  Determination  of 
the  characteristic  functions  of  the  transfer  function  for  the 
updated  system,  (4)  Representation  of  the  transfer  function  for 
the  updated  system  as  a  series,  and  (5)  Laplace  inversion  of  the 
transfer  function. 
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This  report  was  prepared  by  the  AEDAR  Corporation,  under  the  Air 
Force  Office  of  Scientific  Research  contract  number  F49620-88-C- 
0060.  The  work  was  performed  by  Dr.  James  A.  Fabunmi  during  the 
period  May  1,  1988  and  July  31,  1989. 

The  author  wishes  to  thank  Dr.  Anthony  K.  Amos  for  his  support  of 
this  research. 
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I.  OBJECTIVES  OF  THE  RESEARCH 


The  objectives  of  the  research  effort  were: 

(1)  To  develop  MACSYMA-based  routines  for  analyzing  the  Green's 
functions  of  structural  systems  consisting  of  uniform  beams 
supported  at  discrete  locations  by  linear  springs. 

( 2 )  To  extend  these  techniques  to  include  attachments  of  linear 
spring-mass-damper  substructures  to  uniform  beams. 


II.  STATUS  OF  THE  RESEARCH  EFFORT 

After  exploring  different  methods  available  for  the  computer- 
algebraic  representation  of  the  Green's  function  for  flexible 
structures,  a  formulation  was  obtained  which  allowed  the 
automatic  update  of  the  characteristic  parameters  and 
characteristic  functions  by  means  of  which  the  new  Green's 
functions  of  the  modified  structure  can  be  represented.  First  the 
derivations  are  accomplished  on  the  basis  of  the  Laplace 
transform  of  the  Green's  function,  otherwise  known  as  the 
transfer  function,  then  the  inverse  Laplace  transform  of  the 
transfer  function  is  taken  in  order  to  yield  the  Green's 
function.  A  truncated  modal  series  is  used  to  represent  the 
transfer  function. 

The  characteristic  parameters  of  the  updated  transfer  function 
were  obtained  as  the  generally  complex  roots  of  a  characteristic 
polynomial.  Using  the  baseline  parameters  as  a  starting 
approximation,  an  iterative  Newton-Raphson  scheme  was  used  to 
search  for  the  roots  of  the  polynomial.  The  issues  of  convergence 
of  the  truncated  series  and  that  of  the  precision  of  the 
identified  roots  were  investigated.  It  was  found  that  the 
computer  routine  for  performing  these  derivations  can  be 
structured  in  such  a  way  that  any  desired  precision  can  be 
achieved,  at  the  expense  of  retaining  many  terms  in  the  truncated 
series  and  using  multiple  stages  in  the  iterative  search.  Because 
of  limitations  in  the  existing  computational  facility  for 
performing  these  research  tasks,  further  investigation  of  these 
issues  have  been  deferred  to  the  following  year  of  the  effort, 
since  a  dedicated  Symbolics  Workstation  has  been  requested  for 
this  purpose. 

Using  a  general  representation  of  the  Green's  function  of  elastic 
structures,  the  process  of  deriving  the  Green's  function  of  such 
interconnected  structures  consisted  of  the  following  stages:  (1) 
Solution  of  the  integral  equation  which  determines  the  transfer 
function  of  the  updated  system,  (2)  Determination  of  the 
characteristic  parameters  of  the  transfer  Function  for  the 
updated  .system,  (i)  Determination  of  the  characteristic  runctions 
of  the  transfer  function  for  the  updated  system,  (4) 
Representation  of  the  transfer  function  for  the  updated  system  as 
a  series,  and  (5)  Laplace  inversion  of  the  transfer  function. 
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III.  LIST  OF  PUBLICATIONS 


Fabunmi,  J.A.,  "Analysis  of  Modes  and  Frequencies  of  Modified 
Structures  Using  Computer  Algebra",  Proceedings  of  the 
International  Conference  on  Noise  and  Vibration  '89,  Nayang 
Technological  Institute,  Singapore,  August  16-18,  1989. 

(See  Appendix  IV). 


IV.  PROFESSIONAL  PERSONNEL  ASSOCIATED  WITH  RESEARCH  EFirORT 
Dr.  James  A.  Fabunmi. 


V.  INTERACTIONS 

(1)  Clark  Atlanta  University:  Joint  proposal  effort  with  Kwabena 
B.  Bota  -  AFOSR  Proposal  #89-NA-042  entitled  "Analysis  of  Active 
Controller  Effects  on  Flexible  Structures  Using  Computer 
Algebra"  . 


VI.  ADDITIONAL  COMMENTS 

The  formulations  which  have  been  obtained  for  the  interconnection 
of  uniform  beams  to  discrete  substructures  have  turned  out  to  be 
more  general  than  was  originally  anticipated.  Any  one-dimensional 
structure  for  which  the  characteristic  functions  and 
characteristic  parameters  are  available,  can  replace  the  uniform 
beam.  Thus  nonuniform  beams  for  which  a  separate  analysis  and/or 
experimental  tests  have  been  used  to  obtain  this  information,  can 
be  treated.  The  discrete  substructure  is  likewise  not  limited  to 
spring-mass-dampers  only.  Any  attachment  for  which  an  explicit 
transfer  function  is  available  in  the  form  of  a  ratio  of 
polynomials  in  the  Laplace  variable,  can  also  be  handled. 

At  the  moment,  the  theoretical  concepts  for  dealing  with  discrete 
attachments  to  one-dimensional  baseline  structures  appear  to  be 
fully  developed.  The  computer  implementation  of  these  ideas  is 
somewhat  lagging  due  to  the  difficulties  that  have  been 
encountered  in  using  the  MACS YMA  on  the  VAX  via  a  remote 
terminal.  It  is  expected  that  the  dedicated  Symbolics  Workstation 
which  has  been  requested  in  the  follow-on  proposal,  will 
alleviate  these  difficulties  considerably.  Appendices  I  -  III 
are  the  listings  of  the  MACS YMA  codes  which  have  been  developed 
so  far.  Appendix  IV  is  a  copy  of  the  paper  that  was  presented  at 
the  International  Noise  &  Vibration  '89  Conference  in  Singapore. 
This  paper  presents  details  of  the  formulation  of  the 
methodologies  which  have  been  developed  during  this  research 
effort . 
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APPENDIX  I  -  Listing  of  MACSYMA  routine  GF1D.MAC 


comment  ():  =  (  print  l  "  this  module  derives  the  characteristic 
parameter  and  characteristic  function  for  a  specified 
term  in  the  Green's  function  series  for  the  system 
comprising  the  baseline  structure  (batchloaded  with  the 
run)  to  which  the  modification  (batchloaded  with  the 
run)  has  been  attached  at  the  location  specified  in  the 
system  setup  routine,  this  program  assumes  that  the 
baseline  structure  is  one-dimensional  and  the 
modification  is  lumped  parameter,  to  run  the  program, 
type  - 

domoae ( index ) 

where  index  is  the  integer  that  defines  the  desired 
term  in  the  series.  "))$ 

domode ( index ) : =( numer : true , ratprint : false , 

system! ) , idxrange ( index ) , iter ( index ) , print ( "done" ) ) $ 

system (): =block ( print ( "commence  system  setup  routine"), 
againl  ,  bs:read("has  baseline  structure  been  defined? 
*t yes/no )  "  )  , 

if  string ( bs ) ="no"  then  ( load_base 1 ine (), go ( conti  )  ) 
else( 

if  string! bs ) ="yes"  then  go(contl)  else! 
print ( "inappropriate  response" ) , go( againl ) } ) , conti , 
bs: read! "has  sub- structure  attachment  been  defined? 
(yes/no) " ) , 

if  string(bs )="no"  then  ( load_substr ( ) , go ( cont2 ) )  else! 
if  string(bs ) ="yes"  then  go(cont2)  else! 
print! "inappropriate  response" ) ,go( conti ) ) ) ,cont2 , 
print ( "system  setup  complete,  commence  analysis") )$ 

load_baseline( ) : =block( print ( "routine  for  loading  baseline 
structure" ) , 

bsfile : read! "enter  filename  for  baseline  structure 
(extension  .mac  is  assumed )"), load! bsfile ) , 
base_structure ( ) )$ 

load_substr  (  )  :  =  b  1  o  c  k  ( print!  "routine  for  loading 
sub-structure 
attachment" ) , 

ssfile:read( "enter  filename  for  sub-structure 
attachment 

(extension  .mac  is  assumed) "), load! ssfile )  , 
sub_structure { ) , 

bl: read! "enter  attachment  point,  as  fraction  of  length 
of  baseline  structure" ) ,b:bl*el)$ 


idxrange ( argk ) :=block( 

print ( "determine  range  of  indices  needed"), 
startp! argk) , 


3 


P [ 0 ] : pO [ argk ] , val [ 0 ] : p [ 0 ] , 
pz : realpart ( p [ 0 ] ) +%i*imagpart ( p [ 0 ] ) , 
print ("char,  param.  for  baseline  str.  for  this  mode  = 

",  pz), 

icount : 0 , loop , icount : icount+1 , 
itl : argk- icount , itr : argk+ icount , 
if  itlcl  then  itl:l, 
for  kk:itl  thru  itr  do ( startp ( kk ) ) , 
ntr : 1 , newton ( 1 ) ,  val [ icount ] : p [ 1 ] , 
i terms :itr-itl+l, 

print("for  left  index  of  ",itl,"  right  index  of  ",itr, 
"i.e.  no.  of  terms  of  " , iterms ) , 
pz : realpart ( p[ 1 ] ) +%i *imagpart ( p [ 1 ] ) , 
print (" single  iteration  Newton  gives  mod.  char,  param. 
as" ,  pz ) , 

err : abs ( val [ icount ] -val [ icount- 1 ] ) /abs ( val [ icount -1 ] ) , 
print ( "dimensionless  difference  =  " ,err) , 
if  err>0.001  then  go(loop)  )$ 

iter  (  argk  ):  =block  (  print  (  "determine  no.  of  Newton 
iterations" ) , 

ntr : 0 , valpol [ 0 ] : poly ( p[ 0 ] ) , bsl : abs ( valpol [ 0 ] ) + . 001 , 
base :min( bsl ,  1 )  , 
loop, ntr : ntr+1 , 

for  n:l  thru  ntr  do(newton(n) ) , 
valpol [ntr] :poly(p[ntr] ) , 
print ( "valpol  =  ", valpol [ntr ]) , 
err : abs ( valpol [ ntr ] ) /base , 
if  err>0.001  then  go(loop)  else 
pi [ index ] : p [ ntr ] ) $ 

newmode( ) : =block ( pnew: pi [ index] , for  ichkritl  thru  itr 

d  o  (  i  f  pnew  =  pO[ichk]  then 
( oldmode ( ichk ) , go ( exit ) ) )  , 
wc 1 : e v ( wc , p : pnew ) , 

wbpl : sum(phi(b,ks ) **2/ (pnew**2-p0 [ks  ]  **2 )  , 

ks , itl , itr ) , if  ev(wbpl)=0  then 

(brigid( ) , go (exit) ) , 

dnomial : expand! denwc*suml ) , 

hk : ev ( rempoly /dnomial , p : pnew ) , 

if  hk=0  then  (mtplrt( ) ,go(exit) ) , 

for  iold: itl  thru  itr  do( 

top: 2*pnew*phi(b,iold) **2*wcl , 

btm: ( pnew**2-p0 [ iold] **2 ) **2*wbpl*hk, 

alpha2[ index, iold] :top/btm) , 

alpha [ index, iold ] : sqrt ( alpha2 [ index, iold ] ) , 

output ( ) , exit ) $ 

output  (  )  :  =(  print (  "characteristic  parameter 
p ( " , index pnew ) , 

print ("index  of  lowest  term  used  in  series  =",itl), 
print ("index  of  highest  term  used  in  series  =",itr). 
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print ( "number  of  Newton  iterations  used  =",ntr), 
for  imode:itl  tlr— .  itr  do( 

print ( "contribution  of  mode  ",imode,"to  new  mode  =" , 
alpha[ index, imode ] ) ) ) $ 

mtplrt  ():•-(  print  (  "characteristic  polynomial  has  multiple 
roots  for  mode  ",  index,"  this  case  is  not 
yet  considered-  end  of  analysis."))? 

b r i g id ():-(  print (" attachment  point  for  substructure  is 
rigid , 

dynamics  of  baseline  structure  is  unchanged- 
end  of  analysis."))? 

oidmoae ( oldindx ):=( print ( "parameters  of  mode  ", index,"  of 
modified  structure  coincide  with  those  of 
mode  ", oldindx,"  of  baseline  structure- 
end  of  analysis."))? 

poly(pval) :=( prl : product ( p**2-p0 [ kpj * *2 , kp , itl , itr ) , 
suml : sum ( phi (b,ks)**2*prl/(p**2-p0[ks]**2), 
ks , itl , itr ) , 

*  pnomial : expand( denwc*prl-numwc*suml ) , 
nhi : hipow ( pnomial , p ) , 

fnr  nt:0  thru  nhi  do ( a [ nt ]: ratcoef ( pnomial , p , nt )) , 

b[ nhi ] : a [ nhi ] , 

for  nt:nhi-l  step  -1  thru  0 

do ( b [ nt ] :a[nt]+pval*b[nt+l ] )  , 

rempoly : sum( b[ nt ] *p** ( nt-1 ) , nt , 1 , nhi ) , 

b  [  0  ] )  $ 


newton ( n ): =block ( print ( "doing  the  ",n,"th  update  of  a  ",ntr, 
"-round  Newton  iteration"), 
switch: 0 , nn: 0 ,pml : p[ n-1 ] , loop , 
f fprime( n) , 

pp : pml-f [ 0 ] /f [ 1 ] , errp : abs ( pp-pml ) /abs ( pml ) , 

switch: 1 , pml : pp, nn: nn+1 , 

if  errp>.001  then  go (loop), 

print ( nn, "internal  newton  loops"), 

p[n] :pp) $ 

f  fprime  (  n  )  :  =  (  if  n=l  then  (if  switchd  then  routel(  )  else 
route2 ( ) ) 
else  route2 ( ) ) $ 

routel ( ) : =( 

ell:argk,phi(b,ell)  , 
expO : - ( 1/ntr ) *wc*phi**2 , 
f [0] : ev( expO , p : p[ 0 ] )  , 
sumj ( ell , p [ 0 ] ) , 

expl : 2*p[ 0 ] - ( 1/ntr ) *dwc* phi **2+exp0* sumj , 
f [1] : ev{ expl , p: p[ 0 ] ) )$ 
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sumj ( arg] , argp) : =( sumj : 0 , for  j:itl  thru  itr  do( 
if  j=argl  then  term:0  else 
term: 2* argp/ ( argp** 2-pO [ j ] **2 )  , 
sumj :sumj+term) )  $ 

route2 ( ) : = ( 

wbp  :  0  ,  sunk  j  :  0  ,  sunk  :  0  , 

for  iax:itl  thru  itr  do( phi ( b , idx ) , 

df p :  ( pml **  2-pO [ idx ] *  *  2  )  , 

termk: 2*pml/dfp, 

sumk : sumk+termk , 

wbp : wbp+phi *  *  2 /df p , 

sumj ( idx , pml ) , 

surnkj : sumk j +phi *  *  2 *sum j /df p ) , 
expO:l-(n/ntr) *wc*wbp, 
f  [  0 ]  :  ev ( expO , p : pml ) , 

expl : sumk- ( n/ntr ) *dwc*wbp- ( n/ntr ) *wc*sumkj , 
f [ 1 ] : ev ( expl , p : pml ) ) $ 


» 
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APPENDIX  II  -  Listing  of  MACSYMA  routine  FREE_3EAM.MAC 


comment ():=( print ( "this  module  contains  the  functions  needed 
to  specify  the  basis  functions  and  characteristic 
parameters  of  a  freely  suspended  uniform  Euler- 
Bernouili  beam.  The  beam  constants  are  - 
length,  el,  and  EI/mA,  which  are  defined  at  set¬ 
up  time.  The  coordinate  variable,  x,  and  the  mode 
index,  argk,  are  specified  as  calling  arguments  of 

the  functions  phi (x, argk)  and  startp ( argk ) " ) ) $ 

base_structure ( ) :=( comment ( ) , 

el : read ( "enter  beam  length"), 
eima : read ( "enter  value  for  EI/mA"), 
eima2 : sqrt ( eima ) , 
eima 4 : sqrt ( eima 2 ) , 

print ("baseline  structure  is  a  freely  suspended 
Eu le r-Bernoulli  beam,  length  =  "  ,  el, "  EI/mA  = 

" , eima ) ) $ 


phi ( x , argk ) : = ( mu ( argk ) , 

ml:mu*el,mx:mu*x,shx:sinh(mx)  ,  shl :  sinh  ( ml )  , 

sx:sin(mx) ,sl:sin(ml) , chx : cosh ( mx ) , 

chi : cosh (ml) ,cx:cos(mx) , cl : cos ( ml ) , 

snmsl : shl-sl , chmcl : chi -cl , 

shpsl : shl+sl, chpcl : chl+cl , 

shpsx : shx+sx , chpcx : chx+cx , 

phi : (2/ eima 4 ) * ( shmsl* chpcx -chmcl* shpsx ) / 

( shmsl*chpcl-chmcl*shpsl ) ,phi)$ 

startp ( argk ) : =(mu( argk ) ,  val : eima2*mu**2 , 
pO [ argk ] : %i*val ) $ 

mu ( argk ) : =( 

k2  :  %pi * ( argk+1/2 ) , cs : cosh ( k2 ) , 
mu:  (k2*cs-(-l)**argk)/'cs,mu)$ 
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AFFENDIX  ill 


.sting  cf  MACSYMA  routine  MLTHAS . MAC 


lament (): ~ ( print ( "this  module  derives  the  transfer  function  of 
a  local  modification  consisting  of  a  series  of  lumped 
mass/ spring/ dampers  which  may  or  may  not  be  attached  to 
ground,  depending  on  the  values  assigned  to  the  components 
(see  documentation ) .  the  output  of  this  module  are: 

wc  =wc(p)  -  the  impedance  function  of  the  substructure 
as  seen  at  the  point  of  attachment  (driving  point 
impedance ) , 

nunv:  -  the  numerator  cf  the  rational  form  of  wc , 
denwc  -  the  denominator  of  the  rational  form  cf  w c,  and 
awe  -  the  derivative  of  wc  with  respect  to  p  "))S 

sub_s true cure (  ) :  =  ( comment (  ) , input (  ) , process (  ) , printout !  ) S 

input ():=( nd : read ( "enter  order  of  substructure"), 

print ( "enter  substructure  parameters  as  requested"', 
for  id:J  thru  nd  do( 

m [ id ] :read(mass[id] ) , k[ id ] : read ( spring [ id ] ) , 
c [ id ] : read ( damper [ id ] ) ) ) $ 

process! ) : =block ( ratprint : false , 

ndl:nd-l,xf  ndlcl  then  ( xx [ 2 ] : 0 , go ( cont ) ) , 
xx [ nd ) : ( k [ ndl ] +c [ ndl ] *p ) / ( m[ nd ] *p**  2  + 

( c [ ndl ] +c [ nd ] ) *p+k[ ndl ]+kf  nd J ) , 
if  ndl<2  then  go(cont), 
for  id:ndl  step  -1  thru  2  do( 
idl : id- 1 , 

xx [ id ] : ( k [ idl ] +c [ idl ] *p ) / ( m[ id ] *p**2+ 

(c[idl]+c[id]*( l-xx[ id+1 ] ) ) *p+ 
k [ idl ] +k [ id] * ( l-xx[ id+1 ] ) ) ) , cont , 

wc : ratsimp( - ( k [ 1 } +c [ 1 ] *p ) * ( 1 -xx [ 2 ] )-m[l]*p**2,p) , 
dwe : dif  f ( wc , p ) , 
numwc : num ( wc ) , 
denwc  :  denoin  ( wc  )  )  S 

printout ():=( print ( "transfer  function  wc(p)  =  " ,wc) , 
print ( "numerator  of  wc ( p )  =", numwc), 
print ( "denominator  of  wc(p)  =", denwc), 
print! "set  up  of  substructure  transfer 
complete" ) ) $ 


function 
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The  objective  of  this  algorithm  is  to  derive  the 
coefficients  akn  ,  using  information  obtainable  from  the 
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